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Abstract 

Thermostats are dynamical equations used to model thermodynamic variables such 
as temperature and pressure in molecular simulations. For computationally intensive 
problems such as the simulation of biomolecules, we propose to average over fast mo- 
mentum degrees of freedom and construct thermostat equations in configuration space. 
The equations of motion are deterministic analogues of the Smoluchowski dynamics in 
the method of stochastic differential equations. 

1 Introduction 

One of the most ambitious challenges in mathematical modelling of biological processes is to 
describe dynamics of two major biological events within a cell - DNA molecule replication 
and transcription. In both cases, dynamical properties, especially large amplitude confor- 
mational changes of the double-stranded DNA molecule play a vital part. The principal 
physical feature of biological functioning of biomolecules is that they operate at ambient 
physiological temperature, pressure and solvent conditions. Thus, the surrounding physio- 
logical solvent plays the role of a thermostat, among others. To properly thermostat this 
dynamics, especially on a biological ("slow") time scale, we propose a suitable and effective 
temperature control for both the practice of numerical simulation and the general theory of 
dynamical systems. 

For a recent and comprehensive review of the problem outlined above we refer to 
Current approaches commonly use the Nose-Hoover thermostat method (for review see 
m 010113]). This canonical thermostat method involves integration of both position and 
momentum phase space variables. However, for problems that are related to slow conforma- 
tional changes of biomolecules, integration of fast momentum variables appears superfluous 
from a theoretical point of view (unobservable variables) as well expensive in the sense of 
numerical simulation. 

In this Letter, a novel approach to the problem of slow conformational changes in ther- 
mostatting dynamics is presented. The method is based on an analogy with the derivation of 
the Smoluchowski dynamics (Eq. (O below) from the Langevin stochastic dynamics (Eq. 
below) EI ■ The configurational deterministic thermostat is constructed so as to effect 
the temperature control via certain dynamics of the relaxation rate variable and also involves 
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a dynamically fluctuated collective force to ensure the ergodicity property. As the result, the 
temperature control is connected to the specific configurational temperature recently intro- 
duced in 1^^] in a different context. The new configurational thermostat can be combined 
with complementary temperature control via a dynamically fluctuated virial function 
that helps to enhance the efficiency of the thermostat and, more importantly because biolog- 
ical molecules are functioning at a constant pressure, implements the pressure control into 
the corresponding dynamics. Moreover, the configurational thermostat admits stimulation 
by a chain of thermostats similar to Nose-Hoover jj^l and a stochastically driven method 
(see below). To test the new configurational thermostat, corresponding simulations of a 
one-dimensional harmonic and the Morse oscillator dynamics is given, providing a stringent 
test of the ergodicity property. 

Remark: Braga and Travis ^31 have recently proposed a thermostat based on the Smolu- 
chowski equation and configurational temperature ideas (but retaining the momentum vari- 
ables), citing an early preprint of this paper |14| . 

2 Preliminaries 

The early and successful attempts to describe the dynamics of a mechanical system being in 
contact with an environment playing a role of the thermostat is based on the concept of the 
stochastic differential equation [S] . The case of one-dimensional motion of a particle mass 
m in a potential field V{q) provides the characteristic example, 

mq^p, p= -VV{q)-jp + V2Df{t), (1) 

where the friction coefficient 7 and the intensity D of the external random force f{t) are 
connected by the relation, D — m^/k-QT ; here f{t) is the generalized Gaussian stochastic 
process, "white noise", with characteristic cumulants (/(t)) = 0, {f{t)f{t')) — 5{t ~ t'). 
The equilibrium solution of the corresponding Fokker-Planck equation j6j is known, p^o oc 
exp [- (pV2™ -I- F((Z)) /(fcsT)] . 

Langevin's equation is a prototype of the Nose-Hoover deterministic dynamics (after 
generalisation of the last with the term ^ as in Eq. {T)) it is especially evident) . But while 
the Langevin dynamics generate all the sample trajectories and the corresponding measure, 
the Nose-Hoover dynamics produce a single sample trajectory with the correct canonical 
ensemble statistics when ergodicity holds. 

Often details of the dynamics of a system on short time scales are not needed for a 
dynamical description of the observable variables. However, then the Smoluchowski limit of 

o, 

mq = -TWV{q) + ^/M^f{t), (2) 

where r = 7"^, is an appropriate formulation El- The position variable only is 
involved in this equation. Formally, it is supposed that the momentum variable relaxes 
to the local equilibrium state. The corresponding Smoluchowski equation has the Boltz- 
mann distribution as the equilibrium solution 0. Eq. Q, without the random perturba- 
tion, appears as a dissipative dynamics with V{q) playing a role of the Lyapunov function, 
V = —T{yV{q)) < 0. Thus, the full dynamics is a superposition of relaxation to a mini- 
mum of the potential and random perturbations that occasionally expel the system outside 
the vicinity of the minimum. This process equilibrates the system. 
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3 Configurational thermostat 



It is reasonable, in the spirit of deterministic thermostat methods, to conjecture that it is 
possible to use the relaxation time r for thermostatting configurational degrees of freedom 
when momentum variables are still relaxed in their local equilibrium state. Of course, in 
this case the sign of r is not fixed and V loses its meaning as a Lyapunov function. In a 
sense, it means that time can go back as well as forward. 

First, consider the simple dynamical equations for Af degrees of freedom, 

mq = -TVV{q), (3) 

where r is a constant, but by analogy with Nose-Hoover will be endowed with its own 
equation of motion below. By the change of variables, x — ^/mq, it is possible to exclude 
all masses from the formulae in what follows, but we prefer to save the physical notation. 
Requiring that the corresponding Liouville equation has the Boltzmann distribution, p^o oc 
exp [—V{q)/ {kBT)], as a steady state solution, we arrive at the condition that involves the 
temperature in the dynamics, 



E- 



= 0, (4) 



where ^ denotes summation over all particles of the system. After time averaging, Eq. Q 
appears as the definition of the recently introduced so-called configurational temperature 
llOj . Currently it is used in molecular dynamics simulations jl5| . In a more general 
context, in the case of a presupposed anisotropy in the system, let us assume that r in Q 
is a matrix, r ^ F. Then the dynamics take the form, mq = —T'W{q), and the condition 
that involves the temperature in the dynamics is 
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Note that this condition involves the presupposed time scales in the dynamical temperature 
control. Conventional Nose-Hoover methods do not allow such a generalisation. In what 
follows we also consider r as a scalar. On the other hand, it is useful to keep in mind the 
possibility of the generalisation. 

We now attempt to generate statistics as in the Nose-Hoover scheme by making r an 
independent variable in (O . It is easily seen that this is too simple. At an equilibrium point 
'VV = 0, that is, all forces are zero, the evolution comes to a halt and no longer fluctuates, 
irrespective of the time dependence of r. For initial conditions with nonzero forces © after 
a (positive or negative) change of time variable, it is a gradient flow as defined in 1 and it 
is easy to show that all trajectories move along paths in q with equilibrium points at either 
end. In short, the system is not ergodic. Note also that l@J is singular when 'VV = 0. 

The way to overcome the difficulty is suggested by the Smoluchowski stochastic equation. 
In this equation the ergodic motion is ensured by the random forcing. Hence, we need to 
add a deterministic analogue of the random force term in Let us consider, instead of 
l(2Jl, the dynamical equations, 

mci= -TVV{q) + ^, (5) 
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where ^ are constant vectors, but they will be endowed with their own equation of motion 
below. To detail the nature of vectors ^, three principal cases are possible: (a) All ^ = 
{^i}iLiy where N is number of particle in the system, can be varied independently; (b) All 
ii = i Sire varied identically; (c) There exist a preferred direction, e, where e is a constant 
unit vector of physical {i.e. three dimensional) space, and only one variable, C = ^e, 
is varied. Requiring that the Liouville equation corresponding to |SJ| has the Boltzmann 
distribution as a steady states solution and varying ^ according to cases (a), (b) and (c), 
together with temperature control condition we correspondingly obtain the following 
conditions, 

(a) VV{q) = 0, (b) V-Vy(g) = 0, (c) V -e • Vy(g) = 0. (6) 

^ — ' TO ^ — ' TO 

These conditions do not involve temperature but the thermalized forces acting in the system. 
All of them, as well as atheir combination, are candidates for simulating the deterministic 
analogue of the random force, chosen according to the problem under consideration. For 
example, in respect of the Peyrard-Bishop dynamical model of the DNA molecule Q], case 
(c) appears to be appropriate. The physical sense of the conditions above are the following. 
Case (a): the force acting on a particle in the system equals zero (static equilibrium of 
forces). Case (b): the static equihbrium offerees is not required but the total force acting 
on the system equals zero (stability of the system). Case (c): stability of the system in the 
direction e. 

It is practical to remark that a thermostatting dynamics more general than ((SJ is possible, 

TOq = —T\/V{q) + r]mq + ^, (7) 

where the term yymq is suggested by the virial thermostatting scheme that first appears in 
in the context of the harmonic oscillator. Requiring that the corresponding lO Liouville 
equation has the Boltzmann distribution as a steady state solution we obtain together with 
I0J andEJ the following condition on the virial function, 

MkBT - ^ q • = 0. (8) 

In this case a double temperature control is provided. The rj term is not a mandatory tem- 
perature control for our configurational thermostat. We could consider only the t term by a 
trivial modification of 10 and subsequent equations but consider a more general equations 
of the form The reason is that the virial function involves the pressure in the dynam- 
ics and this is important for biologically oriented models. For the sake of dcfiniteness and 
keeping in mind a future application to the Peyrard-Bishop model of the DNA dynamics we 
here fix case (c). When ^ is endowed with its own equation of motion, it provides a shaking 
of the system along the direction e. 

Now we have three parameters which we group together as a 3- vector ct = (r, ?/, S,)"^ . As 
in lO and © we find a stationary solution of the Liouville equation of the form 

Poo (X exp [-{V{q) + a^Qa/2)/{knT)\ 

where Q is a positive definite real symmetric matrix. It is, of course, not possible to 
justify consistently that parameters ot have a Gaussian distribution at equilibrium. It is 
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just a reasonable assumption that leads to the simplest form of self-consistent thermostat 
dynamics. Intuitively, we can justify the Gaussian character of the parameters at equilibrium 
by the limiting theorems of probability theory. Note that when multiple thermostats are 
present, they are usually assumed to be uncoupled (diagonal Q); mathematically it is not 
required and we do not feel that this is physically necessary, thus we consider that the more 
general case of coupled thermostats may be useful. For comparison with the uncoupled case 
we define Q — diag{QrT, Qr/ri, Qcc) to be the matrix with only the diagonal components of 
Q. Eq. Q becomes 

Qg=\ AAfceT-Eq-Vy ^0, (9) 
V 'Eie-VV J 

which now defines g. 

Let us realize the main conjecture the configurational thermostat scheme and allow the 
components of a to fluctuate so that lO holds only after time averaging; we write 

a = G (10) 

where G is as yet an undetermined vector of functions. Now requiring the same condition 
for the solution of the Liouville equation we find that the only solution is 

G = Q^Qg. (11) 

Thus the only undetermined parameters of our thermostatting scheme are the components 
of the positive definite real symmetric matrix Q, and in the uncoupled case Q = Q we have 
G = g. 

We can now ask whether the addition of new variables 77 and/or 1^ will remove the lack 
of ergodicity implied by the potential flow argument applying to (PJ. A partial answer is 
provided by the Frobenius theorem of differential geometry ^7] , which in our case states that 
an integral surface exists (hence the dynamics are definitely not ergodic) if a vector space 
containing the terms in the equation for q but smaller than the full phase space is closed 
under Lie brackets. For realistic potentials (not the harmonic oscillator) this is very unlikely 
since multiple derivatives of V are almost always linearly independent. If the theorem does 
not apply we are in the same situation as for nonthermostatted nonintegrable many particle 
systems, which are often assumed to be ergodic, at least for practical purposes. 

Since under the transformation t —t, a — > —a the equations of motion and Hl()|l 
are still unchanged, they are time reversible. 

To find a mechanically important integral of motion of system (UJ, we need to add a 
redundant variable. Indeed, consider the balance of the mechanical work along trajectories 
of Eqs. and lini), 

^ -WV ■ dq = d{a^Qa.l2) + - Afrj^ dt. 

To obtain an exact differential equation, it is necessary to set 

AV 



T -Mr]]dt^ de. 

m 
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In that case, the foUowing first integral takes place, 

Is = V (q) + + k^TO. 

Since the origin of the redundant variable 9 is arbitrary, it is always possible for an arbitrary 
fixed trajectory to set Is = 0. This integral of motion is useful as a control parameter in 
numerical simulations. Besides, it clearly relates to the equilibrium distribution poo and 
thus can be considered as a first step to reformulation of Eqs. l(7|)- (|10l) in terms of a free 
energy functional. Recall that the corresponding Nose-Hoover integral of motion |21 El El El 
is given by its Hamiltonian - no Hamiltonian is possible here since the momentum does not 
appear explicitly. 

4 Stimulated configurational thermostats 

Since the configurational dynamics above result in the Gaussian equilibrium fluctuation of 
thermostat variables, the latter admits reinforcing by a chain of equations analogous to the 
Nose-Hoover chain thermostat |12) . The chain method consists in including a subsidiary 
sequence of dynamical variables, {cci}, into a thermostat scheme such that asymptotically, 
in the equilibrium distribution, they are independent Gaussian variables, 

Poo (X exp[-(T/(g) + ]-a^Qa + V Q.cxi) / {kBT)]. 

The corresponding dynamics are not unique. We have obtained a clear method for general- 
izing the chain, but since this does not directly relate to our main topic we do not discuss 
the details. Instead we cite the example of the chain rule that has been used in our test 
simulation, 

mq = -TW{q) -|- ^e, f ^ + nr, n = — ^ [k^T - Qr,_^T'}^^) + n+m, ^ = g^, (12) 

where i = 1, . . . , M, tq = t, tm+i = 0. It is a simple chain of total length M. 

It should also be noted that in spite of its popularity, the effectiveness of the chain 
method for computing non-equilibrium properties has been questioned [TH] . 

It is possible to stimulate the Gaussian fluctuation of the thermostat variables by the 
process of Brownian motion. This scheme has the advantage of ensuring the ergodicity 
property. The stimulation, similar to the chain one, may be done in a few ways, applying 
to one or more of the variables t, 77 and ^. In general we have 

mq = -TW{q) +r]mq + ^e, a = G - Aa + V2Df{t), (13) 

where now A and V2D are positive definite real symmetric matrices, and f{t) is a vector 
of independent white noise components. The Liouville equation corresponding to H13|l . 
averaged over all the realizations of f{t), has the form of the Fokker-Planck equation and 
the Boltzmann distribution as a steady state solution only if 

fceTA = DQ. 
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Figure 1: Probability distributions of position variable (shown on background of exact ana- 
lytical distribution) of the harmonic (1) and Morse (2) oscillator. Probability densities are 
calculated as normalized sojourn distributions. Correspondingly, thermostat is: (al)-(a2) 
non-stimulated and uncorrelated (Eq. Q, t — gr,ri = 0,£_ ~ g^); (bl)-(b2) non-stimulated 
and uncorrelated but under double temperature control (Eq. (Q, f — gT,V — g-qA = 
(cl)-(c2) uncorrelated but stimulated by the chain rule (Eq. H12I) . M = 1); (dl)-(d2) stim- 
ulated by stochastic process (Eq. H13() . rj = 0, D = 1). All simulations are performed at 
Qt = Q{ = 1, Qr; = 0.1 for t = lO"* (squares) and t = 10^ (black circles). If squares do not 
shown then distributions are sufficiently near or coincide (al). 

We consider that the most physical case is when the noise is used only for temperature 
control, that is, for t and rj only. We do not establish extreme generality here because our 
main aim is the presentation of the idea of the deterministic configurational thermostat. 
The effectiveness of Eqs. (jlBfl for systems far from equilibrium is not clear. 

5 Test numerical simulations 

The harmonic oscillator is both a simple and an important physical system. At the same 
time, it reveals the ergodicity problem in the canonical ensemble simulation. For this reason, 
it is important to test the Smoluchowski thermostat method capable of generating the 
Boltzmann distribution for a single harmonic oscillator in one dimension, V{q) = /2. 
Then it is reasonable to simulate another one-dimensional system, 'good' from the point 
of view of the Frobenius theorem, and to compare results. We choose for this purpose 
the Morse oscillator, V{q) — Vq{1 — exp(— ag))^ + kg^/2. Simulations are performed using 
global parameters m = 1 and k^T ~ 1, and the Morse potential parameters Vq = 0.25, 
a = 2, k = 0.25. Figure ^ shows the probability distribution of the position variable 
calculated with the four simplest configurational thermostats. Note the effectiveness of the 
double temperature control. The ability of the Smoluchowski thermostat to reproduce the 
correct distribution function, p{q), even with absolute minimum of this thermostat capacity, 
demonstrates its great potential for application. 
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6 Conclusion 



An innovative constant temperature thermostat, the configurational thermostat, exclusively 
involving dynamics of the configurational variables has been introduced. It poses the gen- 
eral problem of the derivation of a thermostatting dynamics for slow dynamical variables 
and outlines at least one way of the solution. For practical purposes, the new thermostat- 
ting scheme can easily be combined with complementary temperature control via dynamical 
fluctuation of the virial function. This combination helps to enhance the efficiency of the 
thermostat temperature control as well as to implement pressure into the dynamics. It is 
also relevant to emphasise the appearance of the dynamically fluctuated forces and corre- 
lated temperature control in the presented thermostatting scheme. We finally remark that 
the proposed method is applicable to thermostatting such macro/meso-scale models as the 
Ginzburg-Landau or reaction-diffusion dynamical equations. 

AS is grateful for support from the University of Bristol and the University of Dundee, 
from CCP5, and from the Royal Society (London). 
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